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ABSTRACT 


A method for constructing upper bound estimates for device 
single event upset (SEU) rates is proposed. A directional Hein- 
rich flux, as a function of direction, must be known. A computer 
code, included in this publication, converts the directional 
Heinrich flux into an "effective flux". The effective flux pro- 
vides a simple way to estimate upper bound SEU rates for devices 
with a known normal incident cross section versus LET curve. 
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1. INTRODUCTION 


Particle sources that can cause single event upset (SEU) in 
integrated circuits include solar events, the anomalous component 
(which is important in some low Earth orbits) , and trapped heavy 
ions around a planet (e.g., Jupiter). For many space projects, 
one or another of these components will dominate the mission 
integrated fluence (which determines the probability of an SEU 
during the mission) and/or peak flux (which impacts the effec- 
tiveness of some error detection and correction schemes) . Thus, 
for many space projects, one or another of the above components 
is the most important SEU related environment. All of these 
components are relatively soft in the sense that mass shielding 
can cause significant attenuation. Galactic cosmic rays are hard 
(only slightly attenuated by mass shielding) but are frequently 
not the environment that causes a space project the most trouble. 
The most trouble often comes from the softer components. Because 
these components are affected by mass shielding, and such shield- 
ing is rarely symmetrically distributed in a spacecraft, the flux 
inside a spacecraft is not isotropic. Therefore, there is a need 
for the ability to estimate SEU rates in a nonisotropic flux. A 
particular calculational method is presented here. The calcula- 
tional method applies to SEU caused by direct ionization. Proton 
induced SEU involving nuclear reaction products requires a sepa- 
rate calculation which is not discussed here. 


The SEU rate for a device can be expressed as 


rate = 


*00 

*1 r 

J 0 

-i. 


27 T 


h(L,0,0) <j (L, 9 , 0) d 0 d(cos9) dL 


( 1 ) 


where h(L,©,0) is the differential (in LET) directional flux 
evaluated at LET L and at the spherical coordinate angles 0 and 0 
(the Z axis is normal to the device) . a is a LET dependent direc- 
tional cross section. 
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The best determination of a typically comes from heavy ion 
tests using particle accelerators. (There are special conditions, 
discussed in section 4, where a can be deduced from a simple 
physical model. The present discussion applies to the more common 
situation where a is known only from heavy ion tests.) If a was 
known in sufficient detail (it is assumed that the flux, h, is 
known) , the rate could be evaluated by numerically evaluating the 
integral in (1). Unfortunately, a is usually only partially 
determined from experimental data because data at large e is not 
meaningful (i.e., does not apply to particles found in space) 
unless the test ions have a lot of range. Machines capable of 
producing such long range ions (e.g., the Bevalac) are expensive 
to use and, therefore, most SEU tests do not use such machines. 
JPL normally tests with lower energy machines and restricts 0 to 
be between 0 and 60°. With 0 so restricted, the "cosine law", 
which states that 

ct(L,©,0) k |cos©| a (L/ |cos0| ,0,0) , (2) 

is usually observed. The large 0 (large enough to violate (2)) 
behavior is usually not measured. 

Because the large 0 behavior of a is usually unknown, it is 
necessary to make worst-case assumptions, i.e., we look for an 
upper bound for the SEU rate. This will be done by assuming a 
particular physical model and adjusting unknown parameters to 
obtain a maximum predicted rate 1 (details are discussed later) . 


l.Some devices exhibit a significant deviation from the cosine 
law for angles as small as 60° or less. For such devices, the 
existing data qualifies as large angle data (the tested angles 
are large enough to violate (2)) and it is possible to justify 
SEU rate estimates that are less pessimistic than the upper bound 
estimates considered here. Instead of adjusting model parameters 
for a maximum predicted rate, they are adjusted for a best fit to 
the measured data. 
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There are circumstances where upper bound estimates can over- 
estimate the actual SEU rate by orders of magnitude. This can 
occur when h is of such a form that the SEU rate can be dominated 
by particles that hit the device at grazing angles (9«90°). In 
this case, the most important angles are the large angles (beyond 
the range such that (2) applies), and typical experimental data 
is inadequate because the test is performed at the least impor- 
tant angles. There are no data pertaining to the most important 
angles, and the data leave SEU rates undetermined by orders of 
magnitude. An estimate that is an upper bound (i.e., large enough 
to compensate for the inadequacy of the data) could be orders of 
magnitude larger than the actual rate. 

Upper bound estimates are occasionally (not always) much higher 
than the actual SEU rate, but such estimates can still be useful 
to a space project. It often happens that a device is hard and/or 
the environment is mild, so that even an excessively pessimistic 
prediction does not bother the project. If this is the case, an 
upper bound estimate is good enough. If not, the project might 
have to obtain and utilize large incident angle cross section 
data so a less pessimistic estimate can be obtained. This was 
done in the past by JPL, but the expense was considerable. 

Because upper bound estimates are useful and more economical 
than accurate methods, they are the final objective of this 
publication. However, as a means to an end, the SEU rate will 
also be expressed as a function of the parameters of a particular 
physical model. This allows the rate to be estimated when the 
parameters happen to be known. Most of the analysis to follow is 
concerned with expressing the rate in terms of model parameters. 
Adjusting the parameters to maximize the rate is discussed in 
Sections 7 and 8. The "maximized rate" can be used to construct 
an "effective flux" which greatly simplifies SEU rate calcula- 
tions (see Sections 7 and 8) . 
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2. THE PHYSICAL MODEL 


We will use the credible (but nonrigorous) assumption that 
there will be enough pessimism in our treatment of the 0 depend- 
ence of a to more than compensate for any possible <p dependence. 
In other words, we assume that an upper bound estimate, derived 
for devices with azimuthal symmetry, will over-estimate an actual 
SEU rate even if the actual device does not have azimuthal symme- 
try. Thus, we confine our attention to devices with azimuthal 
symmetry and cr(L,9) will replace a(L, ©,<£). 

Another assumption is that the device consists of a number of 
charge-collecting volumes and an SEU occurs when one of these 
volumes collects an amount of charge that exceeds some critical 
value. It is also assumed that the charge collected by a volume 
is proportional to ion path length (in the volume) times ion LET 
(this is a traditional assumption [1]). For the reasons given in 
the previous paragraph, we will confine our attention to circu- 
lar-cylindrical charge-collecting volumes. 


3. DISCUSSION OF SINGLE CHARGE-COLLECTING VOLUME 

Devices containing a collection of charge-collecting volumes 
will eventually be considered, but the analysis will begin with a 
device (hypothetical or real) containing a single charge-col- 
lecting volume of uniform thickness. The normal incident cross 
section versus LET curve for such a device is a step function of 
height A and with a threshold LET L t ^. The charge-collecting 
volume is assumed to be a circular cylinder of normal incident 
area A and a threshold LET L t ^ is associated with it. The only 
parameter that cannot be determined from the normal incident 
cross section data is the charge-collecting volume thickness T. 
An upper bound SEU rate estimate is obtained by maximizing the 
predicted rate in T with A and L t ^ held fixed. 
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It is interesting that the "cosine law" becomes exact (i.e., 
(2) becomes exact for all 9) in the limit as T->0 and many people 
erroneously think that this limit gives the maximum SEU rate. So- 
called "proofs" of this assertion are flawed because they fail to 
include particle trajectories that enter the volume through a 
side instead of the upper surface. In reality, a plot of rate 
versus T will show a maximum at some nonzero T. An extreme case 
that illustrates this occurs when L th is so small compared to the 
LET of all particles in the environment that virtually every 
particle that hits the volume will cause an upset. In this case, 
the SEU rate will increase as T increases because more particles 
hit the volume when T is larger. Thus the upper bound SEU rate is 
not obtained from the limiting case as T-+0 1 . 

Determination of the maximum (in T) rate requires that we be 
able to estimate the rate for arbitrary T. This is a nontrivial 
calculation and a simplifying approximation is helpful. For a 
given trajectory direction, we can place an imaginary box— shaped 
volume around the charge-collecting volume which is oriented 
relative to the trajectory direction as shown in Figure 1. Any 
trajectory that hits the charge-collecting volume will also hit 
the imaginary box, and have a path length in the box that is at 
least as large as that in the original volume. The path length 
required for an SEU is the same for the two volumes (both volumes 
are assigned the same normal incident threshold LET and have the 
same thickness) so any particle that causes SEU in the original 
volume will also do so in the imaginary box. Therefore, an upper 
bound for a(L,9) can be estimated by replacing the original 
volume with the imaginary box. The box is conveniently oriented 
so that the mathematical analysis (see Section 5) is essentially 
two dimensional. Note that the box has a different orientation 
for different values of the trajectory coordinate <p. For this 


l.The limiting (T->0) rate often is very nearly equal to the 
maximum (in T) rate. But this condition, when it holds, is a 
property of h and L th and does not apply to all possible h and 
L th (as the above counter-example has shown) . 
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(A) (B) 

TRAJECTORY 



Figure 1. Two views of the charge-collecting volume (solid 
curves) and an imaginary box-shaped volume (dotted curves) . (A) 
is a view from a side with the Z axis and particle trajectory in 
the plane of the paper. (B) is a view from above. 
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reason, the box will be called the "rotating box". 


4. UTILIZATION OF DEVICE FEATURE GEOMETRY 

Until now, it was assumed that the only available device data 
is that obtained from heavy ion tests. Occasionally, the geometry 
of various structures (e.g., depletion regions) is known. It is 
tempting to utilize such information, when it is available, and 
this section discusses when and how such information should be 
utilized. Throughout this discussion, it is assumed that the 
charge-collecting volume is a depletion region (DR) . 

If the DR thickness happens to be known, it is tempting to set 
T equal to this thickness (some people also add a "funnel length" 
but this is a mistake and is discussed again later) . If this is 
valid, it is better than selecting T for a maximum rate because 
it will produce an accurate rate estimate instead of an upper 
bound estimate. However, it is often not valid to set T equal to 
the DR thickness, because the physical model (which assumes 
collected charge to be proportional to LET times path length in a 
geometrically well-defined volume) is totally wrong when diffu- 
sion plays a role. 

It is nontrivial to determine when diffusion plays a role (even 
if device time constants are known) because diffusion is both 
fast and slow. If an ion track is close to a DR, there is ini- 
tially a very large carrier gradient between the high density 
track and the sink-like DR boundary. This large gradient produces 
a very strong diffusion current (fast charge collection) which is 
comparable to the drift current during the recovery stage of 
funneling (virtually all charge collection occurs during the 
recovery stage of funneling [2]). As the track dissipates, charge 
collection via diffusion slows down, but a significant charge can 
be collected quickly. Simulations, using a cylindrical coordinate 
version of PISCES, have shown that an ion track can obliterate a 
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DR that is Sjim away in 0.3ns. Closer tracks can do the same in a 
shorter time. 

Although the physical model is wrong when diffusion plays a 
role, an approximate cosine law (for a limited angle range) is 
compatible with diffusion and the physical model can be used in 
a curve fitting sense. When applied to diffusion, the charge- 
collecting volume dimensions that produce the correct SEU rate 
prediction are not the dimensions of any actual physical struc- 
tures. Therefore, dimensions should be determined from heavy-ion 
cross section data, if large angle data are available, and not 
from device feature geometry. If large angle data are not avail- 
able, the appropriate dimensions are unknown; hence the need for 
upper bound estimates. 

If it is somehow known that diffusion can be neglected for a 
particular device, it should be acceptable to set T equal to the 
DR thickness. Many people erroneously think that a "funnel 
length" should be added to the DR thickness. A recent funneling 
analysis [2] concludes that, in the absence of diffusion, total 
collected charge from an ion track can be calculated from the 
physical model (which states that collected charge is proportion- 
al to path length times LET) by extending each DR linear dimen- 
sion (including the lateral dimensions) by a certain percentage. 
The DR can be replaced with an "extended DR" for the purpose of 
calculating collected charge. This reinforces the concept of 
adding a "funnel length" to the dimensions, but the concept cor- 
rectly calculates only charge collected from a given track inter- 
secting the DR. The concept of an extended DR does not correctly 
calculate the cross section. Trajectories intersecting only the 
extended portions do not result in a collected charge (in the 
absence of diffusion) . The funnel length extension of the DR does 
not contribute to the cross section. Hence an extended DR cannot 
play the role of the charge-collecting volume if our simple 
physical model is to be used. 
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It was concluded above that an extended DR should not be used 
to determine the charge-collecting volume dimensions. However, 
the actual DR (no funnel lengths added) can be used in the ab- 
sence of diffusion. The reason is that the total collected charge 
is proportional to the charge liberated in the DR [2]. If total 
collected charge has a threshold or critical value which results 
in SEU, then so does the charge liberated in the actual DR. The 
latter charge can be calculated by applying the simple physical 
model to the actual DR (no funnel lengths are added) . 


In summary, if the importance of diffusion is unknown, it is 
best not to utilize geometric data because we could be putting 
the "right information into the wrong model". A charge-collecting 
volume is a mathematical construct whose dimensions should be 
determined from large angle heavy ion cross section data. In the 
absence of such data, the appropriate dimensions are unknown; 
hence the need for upper bound estimates. If it is known that 
diffusion can be neglected, the appropriate dimensions for the 
charge— col lecting volume are the actual DR dimensions (funnel 
lengths are not added) . Note that the threshold value of charge 
liberated in the DR, rather than the threshold value of collected 
charge, is used to determine the threshold LET (but this is 
irrelevant if threshold LET is known from heavy ion tests) . 


5 . MATHEMATICAL ANALYSIS FOR SINGLE CHARGE-COLLE CTING VOLUME 1 

The objective is to calculate the SEU rate, as a function of A 
(the normal incident area of the original volume) , T, and L th . R 
is the radius of the original volume and is related to A by 
r= (A/tt) . An upper bound estimate will be obtained by replacing 
the original volume with the rotating box (of width 2R) discussed 


1. Readers that are not interested in the mathematical details can 
skip this section. 
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in Section 3. Throughout this section, the cross section cr(L,e) 
refers to the rotating box rather than the original volume. 

We begin the analysis by assuming that 0 satisfies O<0<9O°, 
which represents particles moving in a downward direction, a is 
expressed as the sum of four partial cross sections 

ct(L, 0) = a-^L,©) + ct 2 (L,0) + a 3 (L,0) + a 4 (L,0) 

where is the cross section associated with the i th 
particle trajectories and the groups are defined 
0e (0 , tt/2 ) ) : 

group 1 = Trajectories that enter the box through the upper 
surface and exit through the lower surface. 

group 2 = Trajectories that enter the box through the upper 
surface and exit through a side. 

group 3 = Trajectories that enter the box through a side and 
exit through the lower surface. 

group 4 = Trajectories that enter the box through a side and 
exit through a side. 

The analysis is slightly different depending on whether 
©€ (0,arctan(2R/T) ) or 0€ (arctan(2R/T) ,tt/ 2) and the two cases are 
considered separately. Group 1 trajectories are illustrated for 
the first case in Figure 2. These trajectories all have the same 
path length, 1 (see Figure 2), and SEU will occur if 1 and parti- 
cle LET, L, are related by 


group of 
by (for 


1 L > L th T 


or 
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Figure 2. Illustration of group 1 trajectories when 
O<0<arctan(2R/T) . 




L > L t ^ cos© . 


(3) 


These trajectories can enter the box through any point in the 
upper surface (of area 4R 2 ) except the section of width 1' (and 
area 2R1') shown in the figure. Hence, the trajectories can enter 
through any point of a surface of area 

4R 2 - 2R 1' = 2R ( 2R - T tan©). 

The cross section (when (3) is satisfied) is the projection of 
this area in the direction of the trajectory. This gives: 


a 1 (L,©) 


If 0 < © < arctan(2R/T) then 
2R (2R - T tan©) cos© if cos© < L/L^ 
0 otherwise. 


If arctan (2R/T) <Q<n/2 , there are no group 1 trajectories and we 
get: 


If arctan(2R/T) < © < rr/2 then 
crjJL,©) = 0. 

Figure 3 illustrates group 2 trajectories for the case 
0<©<arctan ( 2R/T) . SEUs are caused by trajectories below the one 
with length, 1 (see Figure 3) , satisfying 

1 = L,. h T/L (4) 

and such trajectories intersect a surface section of height 1' 
(and area 2R1') shown in the figure, where 

1 ' = T - 1 cos© = T - (L.^ T/L) cos© . 

The area of this surface section is 
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Figure 3. Illustration of group 2 trajectories when 
0<9<arctan(2R/T) . 



2R T (1 - L th cose/L) . 


The cross section is this area projected in the direction of the 
trajectory (if 1* is positive, otherwise the cross section is 
zero) and we get: 

If 0 < 0 < arctan(2R/T) then 

2R T (1 - L.^ COS0/L) sin© if cos0 < L/L^ 

a 2 ®) = 

0 otherwise . 

Figure 4 illustrates group 2 trajectories for the case 
arctan ( 2R/T) <0<7r/ 2 . SEUs are caused by the trajectories to the 
left of the one with length 1 (see Figure 4) satisfying (4) and 
such trajectories intersect a surface section of width 1, (and 
area 2R1 ' ) shown in the figure, where 

1' = 2R - 1 sin© = 2R - (L.^ T/L) sin© . 

The area of this surface section is 

2R ( 2R - L th T sine/L) . 

The cross section is this area projected in the direction of the 
trajectory (if 1* is positive, otherwise the cross section is 
zero) and we get: 

If arctan (2R/T) < © < tt/ 2 then 
2R (2R - L th T sine/L) cos© if sin© < 2RL/ (T L^>.) 

°2 (L/ ©) - 

0 otherwise . 

Group 3 trajectories are similar to group 2 trajectories and 
the result is: 
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Figure 4. Illustration of group 2 trajectories when 
arctan (2R/T) <0<rr/2 . 


a 3 (L,0) 


If 0 < © < arctan(2R/T) then 
2R T (1 - cosQ/L) sin© if cos© < L/L t ^ 

0 otherwise . 


a 3 (L,©) = 


If arctan(2R/T) < © < 7T/2 then 
2R ( 2R - L th T sine/L) cos© if sin© < 2RL/ (T L th ) 
0 otherwise . 


The analysis for group 4 trajectories is similar to that for 
group 1 trajectories and the result is: 


If 0 < © < arctan(2R/T) then 
a 4 (L, ©) = 0 


CT 4 (L, ©) 


If arctan(2R/T) < © < rr/2 then 
2R (T - 2R cot©) sin© if sin© < 2RL/ (T L t ^) 
0 otherwise . 


We next consider © in the range 7r/2<0<Tr. For this range (repre- 
senting particles moving upwards) , the groups are redefined by 
interchanging upper and lower surfaces. It is evident that the 
partial cross sections for this range can be obtained from the 
previous expressions by replacing © with n-Q. 

The symbolism can be shortened by using the unit step function 
S Q (the subscript distinguishes it from a Dirac delta function) 
defined by 




1 if x > 0 
0 if x < 0 . 
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The cross section is obtained by adding the partial cross sec- 
tions and the result is: 


If 0 < 9 < arctan(2R/T) then 
ct (L, 0) = f (L, cos©) <S 0 (L/L th - cos0) . (5a) 

If arctan ( 2R/T) < 0 < tt/ 2 then 
a (L, 0) = f (L, cos©) «S 0 (2R L/ (T L th ) - sin0) . (5b) 

If n / 2 < & < n - arctan (2R/T) then 
a (L, 0) = f (L,cos(7r-0) ) 6 Q ( 2R L/ (T L th ) - sin(7T-0)). (5c) 

If u - arctan (2R/T) < Q < n then 
a (L, 0) = f (L,cos(7T-0) ) S 0 (L/L th - cos(tt- 0)) (5d) 


where f is defined by 

f (L,a)=4R 2 a+2RT(l-a 2 ) 1 / 2 -4RT(L th /L)a(l-a 2 ) 1/2 for ae[0,l]. (6) 

For future convenience, it is helpful to define S(L) by 

S (L) = sqrt* ( 1 - 4R 2 L 2 /(T 2 L th 2 )) (V) 

where the * indicates a modified square root function defined by 


sqrt* (a) 


a 1 / 2 if a > 0 
0 if a < 0 . 


( 8 ) 


It is not difficult to show that 

S Q ( 2R L/(T L th ) - sine) = <S o (cos0 - S(L)) for 0 e [0,tt/2) 

S Q ( 2R L/ (T L th ) -sin(7T-0) ) = S Q (cos (tt- 0) -S (L) ) for tt- 0 e [0,rr/2) 
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so that (5b) and (5c) can be written as 


If arctan(2R/T) < e < u/2 then 

ct(L, 9) = f (L, cos©) S Q (cos9 - S (L) ) . (9a) 

If ?r/2 < 0 < n - arctan(2R/T) then 
a(L,9) = f (L,cos(7t-0) ) S o (cos(7T-0) - S(L)) . (9b) 

Having evaluated a, it is necessary to numerically evaluate the 
integral in (1) . The 0 integral can be carried out separately so 
that (1) can be written as 


rate 


2 rr 


00 7 T 

a(L,0) h (L,0) sin© d© dL 

0 Jo 


( 10 ) 


where 


h a (L,e) 


(27T)~ 1 


'2tx 

h(L,0,0) d0 

0 


(ID 


is the 0 average of h. 

The 0 integration in (10) is carried out by partitioning the 
interval [0 ,tt) into the user specified subintervals 

t a i' a i+l) i 1# • • • / N • 

N and a 1 , . . . , a fj+i are specified by the user but they must 

satisfy 0=a^<a2< • • - <a^ + ^=7r. The 0 average flux, h a , is treated as 
a constant (in 0) over each interval and denotes the flux for 

th ai 

the i tn interval, i.e., 


h ai ( L ) = h a ( L ' 0 ) for 0 c ( a i/ a i+i) • (12) 

Another partitioning of [0 ,tt) consists of the intervals 
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r b i' b i+i) 


i 


1, 2, 3, 4 


where 



b ± = 0 

(13a) 


b 2 = arctan(2R/T) 

(13b) 


b 3 = n/2 

(13c) 


b 4 = n - arctan(2R/T) 

(13d) 


b 5 = w . 

( 13e) 

The interval [0,7r) can be expressed as the double union 

4 N 

= i H 1 jS^i.j 

(14) 

where 

^ijj “ [ b i ,b i+l^ ^ t a j' a j+l) * 

(15) 

A more compact expression for C^j can be obtained by defining 

c ^ j — max {aj ,b^} l — 1< •••/ 4 3 1 » •••/ N 

(16a) 

c 2 .^ j = min{aj +1 , b i+1 } i = 1, 4 j = 1, N 

(16b) 

so that C; ^ can be 
± f J 

expressed as 


C i,j = 

t cl i,j' c2 i,j) lf c2 i,j > cl i,j 

empty if c 2 i j < c 1 ^ . 

(17) 
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The C ^ j intervals do not overlap (as seen by (15)) so any inte- 
gral from 0 to ?r can be expressed as 


N 4 
Z Z 
j=l i=l 




where the integral on ^ is defined by 

± f J 

.2 


'If 1 


if J 


if 1 


if j is not empty 


= 0 if C; ^ is empty 
± f J 


'if 1 


and can be written more compactly as 

5 o( c 2 i,j " cl i,j) 


'if 1 


if} 


if J 


(18) 


(19) 


The © integral in (10) can be evaluated using (12), (15), and 

(18) with the result 


V N 4 

a(L,6) h a (L,0) sin© d© = E h a -i (L) Z cr(L,©) sin© d©. 

Jo J-l i=l JCi.-j 


( 20 ) 


To evaluate the integrals on the right side of (20), use (19) and 
note that i=l uses (5a) , i=2 uses (9a) , i=3 uses (9b) , and i=4 
uses (5d) . The result is 


C 


a(L,6) 
If j 


sin© d© = 



* 0 ((L/L th )-z) dz (21a) 



a(L,©) 

2, j 



$ 0 (z-S(L)) dz 


(21b) 
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r D 2 


a (L , 0 ) sin© de = S Q (c 2 3 f j- cl 3 , j ) jH/Z) 5 Q (z-S(L)) dz (21c) 


'3, j 


3 , 3 


r D 2 


<j(L,©) sin© d© = <5 Q (° 2 4 , j" cl 4 , j ^ ( ( L / L th^ ~ z ^ dz ^ 21d ^ 


'4, j 


where the D's are defined by 


4 , 3 


D^. • = cos(c k .j 4 ) for k - 1,2 i- 1,2 

1 , J x t J 


D ^ * = cos 

1 9 J 


j = 1 , . . . ,N (22a) 
( 7 r— j) for k = 1,2 i = 3,4 j = 1 ,...,N . (22b) 


To evaluate the integrals in (21) , define F by 

F (L, z) = 2R 2 z 2 + R T z (1-z 2 ) 1 / 2 + R T arcsin(z) 

+ (4 / 3 ) R T (L th /L) (l-z 2 ) 3/2 if z e [0,1] (23a) 


F(L,z) = 0 


if z < 0 or z > 1 


(23b) 


so that F is an antiderivative (in z) of f. It is helpful to use 
the following identities, which are valid when a,be[ 0 ,l]: 


f (L , z ) S Q ( c-z) dz = 

<S (c-a) [F (L, c) -F (L, a) ] + 5 Q (c-b) [F (L, b) -F (L, c) ] 


(24a) 


'b 

f (L, z) S Q ( z-c) dz = 

J a 

<5 o (a-c) [F (L, c) -F(L, a) ] + 5 Q (b-c) [F (L,b) -F (L, c) ] . (24b) 

Note that if c<0 or c>l, F(L,c) is not defined in (23a) and, in 
order to avoid computer error, it is necessary to extend the 
domain of F. But c<0 or c>l also implies that the F(L,c) terms 
drop out of the above equations so F(L,c) can be defined in any 
convenient way. This was the motivation for (23b) , which extends 


21 



the domain of F in a particular convenient way. 


Applying (24) to (21) gives 



a(L,e) 

1/ j 


sin© d© = 


+ 


i o<° 2 l,j cl l,j> 
*o( c 2 l, j-® 1 !, j) 


« 0 ((L/L th )-D 2 i (j ) 


(F 3 (L)-F 2 1:j (L) ] 
[F 1 lj (L)-F 3 (L)] 


a(L,©) sin© d© = 

C 2 , j 

fi o( c2 2,j“ cl 2 / j) 5 o( D2 2,j- S ( L )) [F 4 (L)-F 2 2fj (L)) 

+ <f o( c2 2 / j~ cl 2,j) 5 o( Dl 2,j- S ( L )) [F X 2 t j ( L ) “F 4 (L) ] 

a(L,©) sin© d© = 

• C 3, j 

{ o( c2 3,j- cl 3,j) 3 o( Dl 3,j- S < L )) 

+ S o( c2 3,j'° 1 3,j ) i o( D2 3,j" S(L)) 

a(L,©) sin© d© = 

C 4 / j 

5 o^ c 4 , j -c 4,j) ( ^ L / L th^ ~ d 1 4 , j ) [F 3 (L) -F 1 4 , j (L) ] 

+ 6 o( c2 4,j- cl 4, j) *o(( L /L th )-D 2 4/j ) [F 2 4>j (L)-F 3 (L)] 

where the F arrays are defined by 


[F 4 (L) -F^ fj (L)] 
C F 2 3 r j (L) -F 4 (L) 3 


(25a) 


(25b) 


(25c) 


(25d) 
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(26a) 


F 1 i( j (L) = FfL.D^-j) 

F 2 i:j (L) = F(L,D 2 i:j ) (26b) 

F 3 (L) = F(L,L/L th ) (26c) 

F 4 (L) = F (L, S (L) ) . (26d) 

Define the array Gj (L) by 

4 f 

(L) = E a (L, 9) sine d0 (27) 

3 i=l C* ^ 

which can be expressed in terms of the F arrays via (25) and 
(27). Combining (10) with (20) and (27) gives 

N foo 

rate = 2n E h a ^ (L) G^ (L) dL . (28) 

j=l Jo 

As an incidental point, we see from (28) that 27rGj (L) is a 
cross section, associated with the j th 9 interval, for upsets for 

particles with LET equal to L. In an isotropic environment, the 

total effective cross section is 

N 

2 n E G-: (L) . 

j = l 

This "effective cross section" is the quantity that the direc- 
tional flux (in an isotropic environment consisting of particles 
all having the same LET L) must be multiplied by to produce the 
SEU rate. The cross section that the omnidirectional flux (direc- 
tional flux times An) must be multiplied by to produce the rate 
is 
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N 

(1/2) Z G i (L) . 

j=l 


To numerically integrate in LET, select LET values L lf 1^ 
for some user specified M where L 1 is to be small enough so that 
the integral from 0 to L^ is negligible (to avoid indeterminate 
forms, it is required that L 1 >0) and 1^ is to be large enough so 
that the integral from Lj^ to the largest possible particle LET is 
negligible. The rate is approximated by 


rate 


2tt 


N M-l 
Z Z 
j=l k=l 


>L k+l 

• L k 


h aj (L) Gj (L) dL . 


It is required that the L's be sufficiently closely spaced so 
that Gj is nearly a constant in each integral. This produces the 
approximation 


J k+1 


h aj (L) Gj (L) dL ~ 

(1/2) [Gj (L k+1 ) + Gj(L k )] 


Lk+l 


h aj (L) dL = 


(l/2)[Gj(L k+1 ) +G j (L k H [H aj (L k ) - H aj (L R+1 ) ] 
where H a j is the integral (in LET) flux. The rate becomes 


N M-l 


rate = ” iii k2l CG :A+1 + G j.kJ t H aj,k - H aJ , k+1 J (29) 


where 


G j ,k _ G j ^ L k^ ( 3 °) 

and H a j k is a user supplied integral (in LET) directional flux 
evaluated at L k and in the direction corresponding to the j*"* 1 8 
interval and averaged in <p. 
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The equations needed for the actual computations are summarized 
below. The user supplies the L's. The b's are calculated from 
(13) and the a's are supplied by the user and are used to calcu- 
late the c's and D's via (16) and (22). The rate is calculated 
from (29) with H a j, k supplied by the user. Gj, k is obtained by 
combining (30), (27), and (25) with the result 


,j ' k = 2 
*o< c l,j 

+ 5 o< c2 l,j 

+ 5 o( c2 2,j 

,2 


+ 5 o< c ‘ 
+ 6 0 {c< 

+ 5 0 (c : 


2,1 
3, j 
3, j 


” cl i,j) 5 o ( ^ L k/ L th^ ° 2 l,j^ “ F lfj» 

~ cl i,j> 5 o<( L k/ L th> " E)1 l,j ) [F *l,j,k _ F 
" « o 0> 2 2fj - s k) l F \ - F 2,j k] 

- c^.j) ^ 0 ( Dl 2,j - s k> t F S,j,k - F k] 

- cl 3,j> 5 o( Dl 3,j - S k) t F k " F 3,j v] 


,kl 

k ] 


C*«» -s) 


3,1 


■< D2 3,j 


- s k> C F 3,j,k - F V 


+ 5 o< c 4, j 
+ { o( c2 4,j 


- F 


- cl 4 , j ) 5 o(( L k/ L th) - D 4,j) £ F ; 

- cl 4,j) 5 o((V L th) - ° 2 4,j) t F %,j,] 


] 


F k^ 


(31) 


where S k =S(L k ) and is evaluated from (7) with the result 
S k = sqrt*(l - 4R 2 L k 2 /(T 2 L th 2 ) ) 


(32) 


and Fl i,j,k =Fl i,j( L k)' F 2 i,j,k =F 2 i,j^ L k)' F 3 k- r 3 ( L k)' F \ F ^ L k) 

and can be evaluated from (26) and (23) with the result 


Fl i,j,k = 2r2 t? 1 !,!) 2 + R T Dl i,j (1 • (D 1.3 

+ R 


) 2 ) 1/2 


T arcsin (D- 1- ^ j ) + (4/3)R T (L^/I^) ^ ^ 

for 1 = 1, 


3/2 
2 


(33a) 


F\ = 2R 2 (L k /L th ) 2 + R T (V L th> (1 “ (V L th> ) 


2x1/2 


+ R T arcsin 


K. L.U' 2 3/2 

( L k/ L th) + ( 4 /3) R T ^th/ 1 ^ (1 " (L k /L th> ) 

if < L+ 


J th 


F \ - 0 


if L k > L th 


(33b) 


F\ = 2R 2 (s k ) 2 + R T s k (1 - (S k ) 2 ) 1 / 2 + R T arcsin(S k ) 


+ (4/3 ) R T (L^/Lfc) (1 - (S k ) 2 ) 3 / 2 • (33c) 
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6. PROGRAM NT FUR 


The analysis in the last section is performed by the computer 
code NIFUR (Nonisotropic Flux Upset Rate) and the FORTRAN source 
code is listed in Appendix A. 

Special compiling considerations may apply if the code is 
modified by changing the DIMENSION statements. This discussion 
uses terminology applicable to IBM PCs using the Microsoft FOR- 
TRAN Compiler with the large memory model. There is a potential 
problem (if modifications are made) due to the large memory 
reserved for the arrays. The arrays will not all fit in the 
default data segment and must be placed elsewhere. As the code is 
written, this placement is automatic because the largest arrays 
exceed 64K and the default procedure is to give them huge ad- 
dresses and place them outside the default data segment. The 
second largest arrays exceed the default data threshold (32K) and 
are automatically placed outside the default data segment. The 
remaining arrays fit in the default data segment. The problem 
occurs if the array sizes are reduced (by editing the DIMENSION 
statements) in such a way that the placement is not automatic but 
still required because the arrays total to more than 64K. In this 
case, we must request the placement when the code is compiled. 
This is done using the Gt option (applicable to the Microsoft 
compiler) . Specifically, the code is compiled with the command 
"FL/Gt NIFUR. FOR" (the Gt option is not needed if the code is not 
modified) . If a different compiler is used, it is up to the user 
to determine the equivalent option applicable to that compiler 
(if such an option is needed) . The compiled code should be tested 
against the example in Section 10 to verify that it is working 
properly. 

The code requires input files which are created the following 
way. The 0 interval [0 ,tt) is divided up into user selected subin- 
tervals and the directional flux is treated as constant (in the 0 
coordinate) on each subinterval. The user supplies a separate 
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table for each 0 interval, which tabulates flux against LET (the 
smallest LET value used in the flux tabulation must be greater 
than zero). An input file HFLUXOO.DAT must be created, which 
identifies the 0 intervals and LET values to be used in the flux 
tables. Another input file, HFLUX01.DAT, must be created, which 
lists the flux as a function of LET for the first 0 interval (LET 
values are not listed in this file because they are in 
HFLUXOO.DAT). Similarly, HFLUX02.DAT contains the flux for the 
second 0 interval, etc. The flux is required to be a directional 
(corresponding to the appropriate 0 interval) integral in LET 
flux which is averaged in the <p coordinate. The input file for- 
mats are described in detail in the source code comment state- 
ments. If there is confusion, the example in Section 10 might 
help. Device data (A, T, and L th ) are entered via prompts and the 
code output (the upper bound SEU rate estimate) is displayed on 
the terminal. 

In the special case of an isotropic flux, SEU rates for box- 
shaped, charge-collecting volumes can be calculated from the 
chord length distribution function for a box [3], and it is 
interesting to compare this exact calculation to NIFUR results. 
Suppose an area A is used in the NIFUR calculation while the 
exact chord length distribution calculation is performed for a 
box with both lateral dimensions equal to A 1 / 2 (L th and T are the 
same for both calculations) . If T is sufficiently small, the 
NIFUR calculated SEU rate will be about 4 /tt or about 1.27 times 
the chord length distribution calculated rate. The reason is 
that, for small T, the rate is proportional to A. Although NIFUR 
is given the area A, which is assigned to a circular-cylindrical 
volume, the area it actually calculates with is the normal inci- 
dent area of the rotating box, which is 4 A/tt . Hence, the two 
calculations disagree by a factor of 1.27 (if T is not small, the 
exact relation between the two calculations is not simple, al- 
though they are roughly equal) . If desired, the NIFUR results can 
be modified by dividing by 1.27. If this is done, we cannot 
rigorously conclude (from the analysis given here) that the 
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modified rate is an upper bound estimate, but it is "probably" an 
upper bound (at least it is nearly equal to an upper bound) and 
it will conform better (for small T) to more exact calculations. 


7. EFFECTIVE FLUX 


In this section, we confine our attention to hypothetical or 
real devices with the property that the normal incident cross 
section versus LET curve can be approximated by a step function. 
The normal incident behavior of such devices is characterized by 
two parameters; a threshold LET, L th , and a cross sectional area, 
A. The "maximized rate" in a given three dimensional environment 
is defined as the rate predicted by NIFUR with thickness, T, 
chosen to maximize the rate. If a device contains a single circu- 
lar-cylindrical, charge-collecting volume, it is evident that the 
maximized rate is an upper bound for the actual SEU rate because 
it was constructed to be such. 

It is, perhaps, less obvious that the maximized rate (calculat- 
ed from L th and A) is also an upper bound for the SEU rate of a 
device containing any number of identical or nonidentical circu- 
lar-cylindrical, charge-collecting volumes, all having the same 
L.^ and with normal incident areas totaling to A. This is true 
for the following reason. Suppose a device contains a collection 
of such volumes, and the thickness of each volume is selected to 
produce a maximum rate for that volume. The sum of the rates for 
each volume is an upper bound for the device SEU rate. It can be 
shown that this choice of thicknesses will result in all volumes 
being geometrically similar. But for geometrically similar vol- 
umes with the same threshold LET, the rates predicted by NIFUR 
are proportional to the normal incident areas, and the sum of the 
rates (an upper bound for the device rate) is the same as the 
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rate for a single volume of normal incident area A 1 (which is the 
maximized rate). Hence, the maximized rate is an upper bound 
whether the device contains one, or many, identical or nonidenti- 
cal circular-cylindrical volumes all having the same threshold 

LET. 

From the above discussion, it is reasonable to expect that the 
maximized rate, for a given A and L th , is an upper bound for the 
SEU rate of any device having the same A and L th . The maximized 
rate is proportional to A and a normalized rate can be obtained 
by dividing the maximized rate by A. This normalized rate, which 
depends (for a fixed environment) only on L th , will be called the 
"effective flux" evaluated at LET=L th . The effective flux evalu- 
ated at LET=L th and multiplied by A is an upper bound for the SEU 
rate for any device with threshold LET L th and cross section A. 

In summary: 

1. Effective flux evaluated at a given LET is constructed by 
running NIFUR with L th equal to the given LET value and 
selecting an arbitrary A. T is varied to produce a maximum 
rate and the effective flux is this maximum rate divided by 

A. 


2. If a tabulation of effective flux versus LET is avail- 
able, an upper bound SEU rate estimate for a device with 
threshold LET L th and cross section A can be obtained by 
evaluating the effective flux at LET=L th and multiplying it 
by A. 


1 . For example, 10 8 cylinders with radius and thickness ec^al to 1 
um produce the same SEU rate as one cylinder with radius and 
thickness equal to 1 cm. The fact that the one . cy^nder is 
nonphysical (the dimensions are much larger than those of any 
physical structures) is irrelevant. 
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8. PROGRAM EFFLUX 


Upper bound SEU rate estimates are trivial if a tabulation of 
effective flux is available. The only calculation that requires 
any work is the construction of the effective flux which is done 
by the EFFLUX code. 

The FORTRAN source code is listed in Appendix B. The same 
compiling considerations applicable to NIFUR apply to this code 
(see Section 6) and the compiled code should be tested against 
the example in Section 10 to verify that it is working properly. 

This program is the same as NIFUR except that the user does not 
specify values for L th , A, or T . Values are automatically as- 
signed as discussed in the last section. The input files are the 
same as those used by NIFUR and the output is an effective flux 
tabulation which will be in the file EFFLUX. OUT. 

EFFLUX tabulates effective flux against the same LET values 
that the original, user specified, flux was tabulated against. 
However, the effective flux is meaningless (and should not be 
used for SEU rate calculations) at the smallest LET values. The 
reason is that a device with threshold LET equal to the smallest 
tabulated LET value can be upset from particles with a smaller 
LET, and such particles are not included in the original flux 
table (their LETs are below the tabulated range) . The effective 
flux is not meaningful unless it is evaluated at a LET that is 
one or more orders of magnitude greater than the smallest LET in 
the user supplied flux tabulation. If a device has a small 
threshold LET, so that the effective flux must be evaluated at a 

small LET, the user supplied flux tabulation must extend to very 
small LET values. 

As discussed in Section 6, it might be desirable to modify the 
effective flux by dividing it by 1.27. It was not rigorously 
shown that rate estimates obtained from such a modified flux are 
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upper bounds for actual rates, but this is a minor modification 
(compared to other sources of error) and will produce better 
conformity with other rate calculation methods. 


9 . DEVICES WITH SMOOTH CROSS SECTION VERSUS LE T CURVES 

Until now it was assumed that a normal incident device cross 
section versus LET curve is a step function. Effective flux can 
also be used to calculate upper bound SEU rates for smooth cross 
section curves. For a smooth curve, the device is modeled as 
consisting of a collection of charge-collecting volumes which do 
not have the same threshold LET. 

Suppose two devices each contain a collection of circular- 
cylindrical, charge-collecting volumes and, at every LET, the 
normal incident device cross section (sum of the step functions 
associated with individual volumes) of the first device is great- 
er than or equal to that of the second device. There need not be 
any relationship between the volumes in the two devices, other 
than the condition that one cross section curve bounds the other. 
It can be shown that the upper bound SEU rate for the first 
device (obtain by using the method of Section 7 for each charge- 
collecting volume and adding rates) is also an upper bound for 
the second device 1 . This implies that if we are given a cross 
section curve for a device containing an unknown collection of 
charge-collecting volumes, an upper bound rate estimate can be 
obtained by calculating the upper bound rate for any convenient 
hypothetical device having a cross section curve that bounds the 
given curve. In particular, if we are given a smooth curve, such 
as the one in Figure 5, we can construct a staircase curve which 


l.This intuitively obvious statement can be proven rigorously, 
but the proof is too long to include here. 
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CROSS SECTION 



Figure 5. An experimental device cross section versus LET curve 
and a staircase curve that bounds it. 
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bounds the given curve, but is otherwise arbitrary. The staircase 
curve represents a hypothetical device which contains one charge- 
collecting volume for each step in the staircase. The normal 
incident area of a given volume is the height of the step associ- 
ated with that volume and the threshold LET of the volume is the 
LET coordinate of the step (see Figure 5) . After finding the area 
and threshold LET of each volume, an upper bound rate is obtained 
for each volume using the method of Section 7 and the rates are 
added to produce an upper bound estimate for the actual device. 
An example is given in the next section to illustrate this algo- 
rithm. 

It is significant that the steps in the staircase curve do not 
have to correspond to physical structures of the actual device. 
Any staircase function that bounds the cross section curve can be 
used, but the lowest (most accurate) upper bound estimates are 
produced by the curves that conform most closely to the actual 

curve . 

As an incidental point, the above algorithm has the mathemati- 
cal interpretation of a numerical integration. The algorithm 
produces an upper bound for the integral 

EF (L) dX (L) 


where EF(L) is the effective flux evaluated at LET-L and X(L) is 
the normal incident device cross section evaluated at LET=L . The 
limiting case of an infinitely fine division (the staircase curve 
has an infinite number of steps and conforms to the actual curve) 
produces the integral exactly. It is not necessary to be aware of 
this in order to perform the calculations, but this information 
is helpful when looking for similarities between the SEU rate 
prediction method discussed here and other methods reported in 
the 1 iter atur e . 
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10. AN EXAMPLE 


The example below estimates the SEU rate for a hypothetical 
device in a hypothetical spacecraft exposed to a hypothetical 
solar flare. 

The mass distribution, relative to the device orientation, in 
this hypothetical spacecraft requires that a different shield 
thickness be used to calculate the flux for each of the angular 
regions shown below: 


Region 

1: 

0 

< 

e 

< 

45°, 

0 

< 

0 

< 

360° 

Region 

2 : 

45° 

< 

e 

< 

75° , 

0 < 

45 

0 

or <p 

Region 

3: 

45° 

< 
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15 °, 
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< 

0 

< 

135° 

Region 

4 : 

45° 

< 

e 

< 

15 ° , 

135° 

< 

0 

< 

225° 

Region 

5: 

45° 

< 

0 

< 

15 °, 

225° 

< 

0 

< 

315° 

Region 

6: 

75° 

< 

0 

< 

90° , 

45° 

< 

0 

< 

135° 

Region 

7: 

75° 

< 

0 

< 

90°, 

0 < 

45 

0 

or <p : 

Region 

8: 

75° 

< 

0 

< 

KD 

O 

o 

225° 

< 

0 

< 

315° 

Region 

9: 

75° 

< 

0 

< 

90°, 

135° 

< 

0 

< 

225° 

Region 

10: 

90° 

< 

0 

< 

180°, 

0 

< 

0 

< 

360° 


where 0 is measured from the device normal and 0 is measured from 
some line lying in the device plane. 

The Heinrich flux behind the shield thicknesses associated with 
each angular region is calculated by a computer code such as 
CREME [1]. A particular code, using a particular (and severe) 
solar flare model and particular shield thicknesses produced the 
flux tables shown in Table l 1 . 


l.This calculation is not discussed in detail here because the 
user must use his own resources (e.g., CREME [1]) to obtain these 
flux tables. The present discussion is concerned with what to do 
with the tables, not how to get them. 
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LET 

. 1000E+0 
. 1259E+0 
. 1585E+0 
. 1995E+0 
. 2512E+0 
. 3 162E+0 
. 3981E+0 
. 5012E+0 
. 6310E+0 
. 7943E+0 
. 1000E+1 
. 1259E+1 
. 1585E+1 
. 1995E+1 
. 2512E+1 
. 3162E+1 
. 3981E+1 
. 5012E+1 
. 6310E+1 
. 7943E+1 
. 1000E+2 
. 1259E+2 
. 1585E+2 
. 1995E+2 
. 2512E+2 
. 3162E+2 
. 3981E+2 
. 5012E+2 
. 6310E+2 
. 7943E+2 
. 1000E+3 


REGION 1 
FLUX 


. 237E+6 
. 161E+6 
. 113E+6 
. 790E+5 
. 567E+5 
. 418E+5 
. 299E+5 
• 2 17E+5 
. 180E+5 
. 152E+5 
. 135E+5 
. 964E+4 
. 660E+4 
. 399E+4 
. 256E+4 
. 179E+4 
. 908E+3 
. 663E+3 
. 351E+3 
. 2 03E+3 
. 121E+3 
. 807E+2 
. 218E+2 
. 110E+2 
. 380E+1 
. 895E-1 
. 323E-1 
. 159E-1 
. 642E-2 
. 240E-2 
. 829E-4 


REGION 2 
FLUX 


. 134E+6 
. 969E+5 
. 708E+5 
. 525E+5 
. 396E+5 
. 300E+5 
. 222E+5 
. 169E+5 
. 142E+5 
. 120E+5 
. 108E+5 
. 755E+4 
. 511E+4 
. 303E+4 
. 192E+4 
. 131E+4 
. 638E+3 
. 464E+3 
. 245E+3 
. 142E+3 
. 850E+2 
. 57 1E+2 
. 156E+2 
. 787E+1 
. 270E+1 
. 559E-1 
. 188E-1 
. 892E-2 
. 352E-2 
. 129E-2 
. 425E-4 


REGION 3 
FLUX 


. 182E+7 
. 117E+7 
. 757E+6 
. 472E+6 
. 303E+6 
. 206E+6 
. 132E+6 
. 813E+5 
. 640E+5 
. 497E+5 
. 403E+5 
. 286E+5 
. 203E+5 
. 136E+5 
. 974E+4 
. 726E+4 
. 472E+4 
. 354E+4 
. 216E+4 
. 130E+4 
. 798E+3 
. 529E+3 
. 145E+3 
. 726E+2 
. 248E+2 
. 354E+0 
. 107E+0 
. 516E-1 
. 208E-1 
.771E-2 
. 285E-3 


REGION 4 
FLUX 


. 325E+6 
. 218E+6 
. 147E+6 
. 100E+6 
. 703E+5 
. 511E+5 
. 357E+5 
. 253E+5 
. 207E+5 
. 172E+5 
. 150E+5 
. 108E+5 
. 741E+4 
. 449E+4 
. 291E+4 
. 204E+4 
. 101E+4 
. 727E+3 
. 385E+3 
. 222E+3 
. 132E+3 
. 880E+2 
. 238E+2 
. 120E+2 
. 4 14E+1 
. 106E+0 
. 395E-1 
. 197E-1 
. 799E-2 
. 299E-2 
. 980E-4 


Table l. Integral directional Heinrich flux versus LET 
angular region. LET in units of MeV-cnr/mg and flux in 
particles/m^-sec-sr . 


REGION 5 
FLUX 


. 224E+6 
. 153E+6 
. 106E+6 
. 750E+5 
. 544E+5 
. 403E+5 
. 290E+5 
. 212E+5 
. 177E+5 
. 149E+5 
. 132E+5 
. 951E+4 
. 650E+4 
. 392E+4 
. 251E+4 
. 175E+4 
. 898E+3 
. 654E+3 
. 346E+3 
. 202E+3 
. 121E+3 
. 806E+2 
. 2 17E+2 
. 110E+2 
. 379E+1 
. 882E-1 
. 3 16E-1 
. 153E-1 
. 625E-2 
. 237E-2 
. 8 19E-4 


for each 
units of 
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REGION 6 
FLUX 


REGION 7 
FLUX 


REGION 8 
FLUX 


REGION 9 
FLUX 


LET 

. 1000E+0 
. 1259E+0 
. 1585E+0 
. 1995E+0 
. 2512E+0 
. 3162E+0 
. 3981E+0 
. 5012E+0 
. 63 10E+0 
.7943E+0 
. 1000E+1 
. 1259E+1 
. 1585E+1 
. 1995E+1 
. 2512E+1 
. 3162E+1 
. 3981E+1 
. 5012E+1 
. 6310E+1 
. 7943E+1 
. lOOOE+2 
. 1259E+2 
. 1585E+2 
. 1995E+2 
. 2512E+2 
. 3 162E+2 
. 3981E+2 
. 5012E+2 
. 63 10E+2 
. 7943E+2 
. 1000E+3 


. 182E+7 
. 117E+7 
. 757E+6 
. 472E+6 
. 303E+6 
. 206E+6 
. 132E+6 
. 813E+5 
. 640E+5 
. 497E+5 
. 403E+5 
. 286E+5 
. 2 03E+5 
. 136E+5 
. 974E+4 
. 726E+4 
. 472E+4 
. 354E+4 
. 216E+4 
. 130E+4 
. 798E+3 
. 529E+3 
. 145E+3 
. 726E+2 
. 248E+2 
. 3 54E+0 
. 107E+0 
. 516E-1 
. 2 08E-1 
. 77 IE-2 
. 285E-3 


. 372E+6 
. 244E+6 
. 164E+6 
. 111E+6 
. 775E+5 
. 561E+5 
. 391E+5 
. 276E+5 
. 225E+5 
• 186E+5 
. 162E+5 
.117E+5 
. 802E+4 
. 487E+4 
. 3 18E+4 
. 2 19E+4 
. 108E+4 
. 778E+3 
. 4 13E+3 
. 236E+3 
. 141E+3 
. 932E+2 
. 252E+2 
. 127E+2 
. 438E+1 
. 118E+0 
•446E-1 
. 220E-1 
. 897E-2 
. 341E-2 
. 123E-3 


. 978E+6 
. 620E+6 
. 400E+6 
. 257E+6 
. 169E+6 
. 117E+6 
. 775E+5 
. 500E+5 
. 399E+5 
. 315E+5 
. 2 61E+5 
. 186E+5 
. 130E+5 
. 828E+4 
. 563E+4 
. 399E+4 
. 228E+4 
. 164E+4 
. 962E+3 
• 579E+3 
. 351E+3 
. 238E+3 
. 635E+2 
. 3 19E+2 
. 109E+2 
. 191E+0 
. 721E-1 
. 353E-1 
. 144E-1 
. 540E-2 
. 199E-3 


. 102E+6 
. 763E+5 
. 572E+5 
. 434E+5 
. 333E+5 
. 255E+5 
. 190E+5 
. 145E+5 
. 122E+5 
. 104E+5 
•931E+4 
. 645E+4 
. 431E+4 
. 254E+4 
. 160E+4 
. 109E+4 
• 542E+3 
. 400E+3 
. 214E+3 
. 124E+3 
. 737E+2 
. 492E+2 
. 132E+2 
. 667E+1 
. 228E+1 
. 401E-1 
. 130E-1 
. 622E-2 
. 249E-2 
. 915E-3 
. 301E-4 


REGIONIO 

FLUX 


. 362E+6 
. 237E+6 
. 160E+6 
. 109E+6 
. 761E+5 
. 549E+5 
. 381E+5 
. 2 68E+5 
• 2 18E+5 
. 181E+5 
. 157E+5 
. 113E+5 
. 777E+4 
. 472E+4 
. 308E+4 
. 213E+4 
. 105E+4 
. 759E+3 
. 402E+3 
. 231E+3 
. 137E+3 
. 908E+2 
. 246E+2 
. 124E+2 
. 427E+1 
. 114E+0 
. 42 9E-1 
. 209E-1 
. 862E-2 
. 329E-2 
. 117E-3 


Table l. (continued) 
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For each 0 interval, the fluxes must be averaged in 0. For each 
0 interval, each 0 interval is the same size so a 0 average is a 
simple numerical average (if the <p intervals were not the same 
size, the fluxes would be weighted according to the sizes of the 
0 intervals in the obvious way) . The Region 1 flux is the <p 
average flux for the 0 interval (0,45°) and the Region 10 flux is 
the 0 average flux for the 0 interval (90°, 180°). The fluxes for 
Regions 2,3,4, and 5 are averaged to produce the 0 average flux 
for the 0 interval (45°, 75°) and the fluxes for Regions 6,7,8, 
and 9 are averaged to produce the 0 average flux for the 0 inter- 
val (75°, 90°) . This averaging is a lot of work if done by hand 
(especially if the 0 intervals are not uniform and a weighted 
average is required) and the user will probably want to write a 
code that will speed this up. The averaged fluxes for our example 
are shown in Table 2. This table is used to produce the input 
files used by EFFLUX as described in the code comment statements. 
The input files are shown in Table 3. 


Before proceeding, this is a good time to test NIFUR. If NIFUR 
is run with the Table 3 files together with A=10 6 , L th =5, and 
T=io~’^ , it should respond with RATE=. 282754 . 


The next step is to run EFFLUX. The output of EFFLUX is the 
effective flux and is shown in Table 4. Because the smallest LET 
in the original flux tabulation is 0.1 MeV-cm 2 /mg, the effective 
flux tabulation should not be trusted for LET less than 1 
MeV-cm 2 /mg (see Section 8). Discarding the lower LET entries and 
dividing the flux by 1.27 (as suggested in Section 8) produces 
the final effective flux table which is shown in Table 5. 


Now that the effective flux has been calculated, the next step 
is to estimate the SEU rate. In this example, the device is 
represented by the following cross section data: 
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0<©<45 ° 

45°<e<75° 

75°<0<9O° 

LET 

FLUX 

FLUX 

FLUX 

. 1000E+0 

. 237E+6 

. 626E+6 

. 818E+6 

. 1259E+0 

. 161E+6 

• 409E+6 

. 528E+6 

. 1585E+0 

. 113E+6 

. 270E+6 

. 345E+6 

. 1995E+0 

. 790E+5 

. 175E+6 

. 221E+6 

. 2512E+0 

. 567E+5 

. 117E+6 

. 146E+6 

. 3162E+0 

. 4 18E+5 

. 819E+5 

. 101E+6 

. 3981E+0 

. 299E+5 

. 547E+5 

. 669E+5 

. 5012E+0 

. 2 17E+5 

. 362E+5 

. 434E+5 

. 6310E+0 

. 180E+5 

. 292E+5 

. 347E+5 

. 7943E+0 

. 152E+5 

. 235E+5 

. 276E+5 

. 1000E+1 

. 135E+5 

. 198E+5 

. 230E+5 

. 1259E+1 

. 964E+4 

. 141E+5 

. 163E+5 

. 1585E+1 

. 660E+4 

. 983E+4 

. 114E+5 

. 1995E+1 

. 399E+4 

. 626E+4 

. 732E+4 

. 2512E+1 

. 256E+4 

. 427E+4 

. 504E+4 

. 3 162E+1 

. 179E+4 

. 309E+4 

. 363E+4 

. 3981E+1 

. 908E+3 

. 182E+4 

-216E+4 

. 5012E+1 

. 663E+3 

. 135E+4 

. 159E+4 

. 63 10E+1 

. 351E+3 

. 784E+3 

. 937E+3 

. 7943E+1 

. 203E+3 

. 467E+3 

. 560E+3 

. 1000E+2 

. 121E+3 

. 284E+3 

. 34 1E+3 

. 1259E+2 

. 807E+2 

. 189E+3 

. 227E+3 

. 1585E+2 

. 218E+2 

. 515E+2 

. 617E+2 

. 1995E+2 

. 110E+2 

. 259E+2 

. 310E+2 

. 2512E+2 

. 380E+1 

. 886E+1 

. 106E+2 

. 3162E+2 

. 895E-1 

. 151E+0 

-176E+0 

. 3981E+2 

. 323E-1 

. 492E-1 

. 592E-1 

. 5012E+2 

. 159E-1 

. 239E-1 

. 288E-1 

. 6310E+2 

. 642E-2 

. 964E-2 

. 117E-1 

. 7943E+2 

. 240E-2 

. 359E-2 

• 436E-2 

. 1000E+3 

. 829E-4 

. 127E-3 

. 159E-3 

Table 2. 0 

average of the 

Table 1 fluxes for the 

vals . The units are the same as in Table 

1 . 


90°<e<180° 

FLUX 


. 362E+6 
. 237E+6 
. 160E+6 
. 109E+6 
. 761E+5 
. 549E+5 
. 381E+5 
. 2 68E+5 
. 2 18E+5 
. 181E+5 
. 157E+5 
. 113E+5 
. 777E+4 
. 472E+4 
. 308E+4 
. 2 13E+4 
. 105E+4 
. 759E+3 
. 402E+3 
. 231E+3 
. 137E+3 
. 908E+2 
. 246E+2 
. 124E+2 
. 427E+1 
. 114E+0 
. 429E-1 
. 209E-1 
. 862E-2 
. 329E-2 
. 117E-3 


four 9 inter- 
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HFLUXOO.DAT 

4 

31 

45.0 

75.0 

90.0 

. 1000E+00 
. 1259E+00 
. 1585E+00 
. 1995E+00 
. 2512E+00 
. 3162E+00 
. 3981E+00 
. 5012E+00 
. 6310E+00 
. 7943E+00 
. 1000E+0 1 
. 1259E+01 
. 1585E+01 
. 1995E+01 
. 2512E+01 
. 3 162E+01 
. 398 1E+01 
. 5012E+01 
. 6310E+01 
. 7943E+01 
. 1000E+02 
. 1259E+02 
. 1585E+02 
. 1995E+02 
. 2512E+02 
. 3 162E+02 
. 3981E+02 
. 5012E+02 
. 63 10E+02 
. 7943E+02 
. 1000E+03 


HFLUX01.DAT 
. 237E+06 
. 161E+06 
. 113E+06 
. 790E+05 
. 567E+05 
. 418E+05 
. 299E+05 
. 2 17E+05 
. 180E+05 
. 152E+05 
. 135E+05 
. 964E+04 
. 660E+04 
. 399E+04 
. 256E+04 
. 179E+04 
. 908E+03 
. 663E+03 
. 351E+03 
. 203E+03 
. 121E+03 
. 807E+02 
. 218E+02 
. 110E+02 
. 380E+01 
. 895E-01 
. 323E-01 
. 159E-01 
. 642E-02 
. 240E-02 
. 829E-04 


H FLUX 02 . DAT 
. 626E+06 
. 409E+06 
. 270E+06 
. 175E+06 
. 117E+06 
. 8 19E+05 
. 547E+05 
. 3 62E+05 
. 292E+05 
. 235E+05 
. 198E+05 
. 141E+05 
. 983E+04 
. 626E+04 
. 427E+04 
. 3 09E+04 
. 182E+04 
. 13 5E+04 
. 784E+03 
. 467E+03 
. 284E+03 
. 189E+03 
. 515E+02 
. 259E+02 
. 886E+01 
. 151E+00 
. 492E-01 
. 239E-01 
. 964E-02 
. 359E-02 
. 127E-03 


HFLUX03.DAT 
. 818E+06 
. 528E+06 
. 345E+06 
. 221E+06 
. 146E+06 
. 101E+06 
. 669E+05 
. 434E+05 
. 347E+05 
. 276E+05 
. 230E+05 
. 163E+05 
. 114E+05 
. 732E+04 
. 504E+04 
. 363E+04 
. 216E+04 
. 159E+04 
. 937E+03 
. 560E+03 
. 341E+03 
. 227E+03 
. 617E+02 
. 310E+02 
. 106E+02 
. 176E+00 
. 592E-01 
. 288E-01 
. 117E-01 
. 436E-02 
. 159E-03 


HFLUX04.DAT 
. 362E+06 
. 237E+06 
. 160E+06 
. 109E+06 
. 761E+05 
. 549E+05 
. 3 8 1E+05 
. 2 68E+05 
. 2 18E+05 
. 181E+05 
. 157E+05 
. 113E+05 
. 777E+04 
. 472E+04 
. 308E+04 
. 2 13E+04 
. 105E+04 
. 759E+03 
. 402E+03 
. 2 3 1E+03 
. 137E+03 
. 908E+02 
. 24 6E+02 
. 124E+02 
. 427E+01 
. 114E+00 
. 429E-01 
. 209E-01 
. 8 62E-02 
. 329E-02 
. 117E-03 


Table 3. Input files. The file names at the top are not part of 
the file contents. 
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EFFECTIVE FLUX 


LET IN UNITS OF MeV-cm**2 /mg 
FLUX IN UNITS OF 1/ (cm**2-day ) 


LET 

FLUX 

0 . 1000E+00 

0 . 5248E+08 

0. 1259E+00 

0 . 4234E+08 

0. 1585E+00 

0 . 3263E+08 

0. 1995E+00 

0 . 2487E+08 

0 . 2512E+00 

0. 1865E+08 

0. 3162E+00 

0. 1387E+08 

0 . 3981E+00 

0 . 1024E+08 

0.5012E+00 

0 . 7514E+07 

0 . 63 10E+00 

0 . 5511E+07 

0 . 7943E+00 

0 . 4078E+07 

0. 1000E+01 

0 . 3064E+07 

0. 1259E+01 

0. 2318E+07 

0. 1585E+01 

0. 1735E+07 

0. 1995E+01 

0 . 1276E+07 

0 . 2512E+01 

0 . 9207E+06 

0. 3162E+01 

0. 6598E+06 

0. 3981E+01 

0 . 4672E+06 

0 . 5012E+01 

0 . 3264E+06 

0 . 63 10E+01 

0 . 2271E+06 

0 . 794 3E+01 

0. 1561E+06 

0. 1000E+02 

0. 1064E+06 

0. 1259E+02 

0 . 7 171E+05 

0. 1585E+02 

0 . 4733E+05 

0. 1995E+02 

0 . 3071E+05 

0 . 2512E+02 

0. 1969E+05 

0 . 3162E+02 

0 . 1249E+05 

0. 3981E+02 

0 . 7847E+04 

0 . 5012E+02 

0 . 4902E+04 

0 . 63 10E+02 

0 . 3058E+04 

0.7943E+02 

0. 1909E+04 

0 . 1000E+03 

0. 1188E+04 


Table 4. Contents of the file EFFLUX. OUT. 
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EFFECTIVE FLUX 


LET IN UNITS OF MeV-cm**2/mg 
FLUX IN UNITS OF 1/ (cm**2-day) 


LET 

0. 1000E+01 
0. 1259E+01 
0. 1585E+01 
0. 1995E+01 
0.2512E+01 
0. 3162E+01 
0 . 3981E+01 
0 . 5012E+01 
0 . 63 10E+01 
0 . 7943E+01 
0. 1000E+02 
0. 1259E+02 
0. 1585E+02 
0. 1995E+02 
0 . 2512E+02 
0. 3162E+02 
0. 3981E+02 
0.5012E+02 
0 . 63 10E+02 
0. 7943E+02 
0. 1000E+03 


FLUX 

0 . 241E+07 
0. 183E+07 
0 . 137E+07 
0. 100E+07 
0 . 725E+06 
0. 520E+06 
0 . 368E+06 
0 . 257E+06 
0. 179E+06 
0 . 123E+06 
0 . 838E+05 
0 . 565E+05 
0 . 373E+05 
0 . 242E+05 
0 . 155E+05 
0 . 983E+04 
0 . 618E+04 
0 . 386E+04 
0.241E+04 
0. 150E+04 
0 . 935E+03 


Table 5. The final effective flux table. 
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Cross section (cm 2 ) LET fMeV-cm 2 /mcM 


Approximately zero 3 

9.2X10 -4 6 

5 . 7xl0~ 3 16 

1 . 6xl0 -2 > 40 


A staircase function is constructed to bound the cross section 
curve as illustrated in Figure 5. The vertices are arbitrary and 
in this example are chosen to correspond to the discrete points 
listed above. Each step in the staircase function defines a cross 
section and threshold LET for a hypothetical charge-collecting 
listed below (note that a given cross section corre- 
sponds to a smaller LET than in the original data, e.g., 9.2xl0 -4 
corresponds to 3 instead of 6) . Also listed are the approximate 
effective fluxes evaluated at the LETs. Upset rate estimates are 
constructed from this data as shown (the obvious units are im- 
plied) . 


Charge-collecting 
volume # 

Cross 

section 

Threshold 

LET 

Effective 

flux 

1 

9 . 2xl0~ 4 

3 

5.7xlO S 

2 

5.7xl0 -3 
-9.2X10" 4 
=4 . 8xl0~ 3 

6 

2.0xl0 S 

3 

1. 6xl0 -2 
-5.7X10 -3 
=1. OxlO -2 

16 

3 . 7xl0 4 

Upper bound SEU rate = 
(9.2xl0“ 4 ) (5.7xl0 5 ) + 

(4.8X10 -3 ) (2 

.OxlO 5 ) + (1. 

OxlO -2 ) (3.7x10 


= 1. 9xl0 3 /day . 
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11. CONCLUSION 


The analysis given here provides a way to standardize SEU rate 
calculations for a given space mission. The proposed method is 
the following: The flight project is required to estimate the 
environments of concern (e.g., typical flux, peak flux, mission 
integrated fluence, etc.) with the effects of mass shielding 
included. Using the computer code EFFLUX, the flight project 
constructs tables of effective flux versus LET. If the environ- 
ment must be treated as nonisotropic, different tables are needed 
for different device orientations relative to the spacecraft. A 
device vendor can then calculate SEU rates by using a simple 
numerical integration to combine device cross section data with 
the effective flux tables. 

The purpose of the effective flux is to simplify SEU rate 
calculations from the vendor's point of view. Most of the work 
required for SEU rate estimates is performed by the flight 
project when constructing effective flux tables. After these 
tables have been constructed, SEU rate estimates are simple. 
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APPENDIX A 


PROGRAM NIFUR 

C Copyright (c) 1991, California Institute of Technology. 

C U.S. Government Sponsorship under NASA Contract NAS7-918 
C is acknowledged. 

REAL LTH , L ( 2 0 0 ) 

CHARACTER* 2 , FLNM1 
CHARACTER* 11 , FLNM2 

DIMENSION A ( 31) ,B(5) , Cl (4, 30) ,C2(4,30) ,D1(4,30) ,D2(4,30) 
DIMENSION S ( 200) , HA (30,200) 

DIMENSION F3 (200) ,F4(200) ,F1(4,30,200) ,F2(4,30,200) 
DIMENSION G ( 30 , 200) 

C Read data. Device data consists of a normal incident cross 
C section AC, a threshold LET LTH, and a thickness T, which 
C are entered via prompts. The files HFLUX00.DAT, HFLUX01.DAT, 
C ..., HFLUX09.DAT, HFLUX10.DAT, HFLUX11.DAT, ... contain 
C environmental data. There is one file (after the "00" file) 

C for each angular bin (discussed below) . The angular inter- 
C val [0,180 degrees) is partitioned into the bins 
C [ A ( 1) , A ( 2 ) ) , [A(2) ,A(3) ) , . . . , [A (N) , A(N+1) ) where A(1)=0 

C and A (N+l) =180 . The discrete values L(l) , ...,L(M) of LET 

C (in MeV-cm**2/mg) are selected by the user. The file 
C HFLUX00.DAT contains the values of N, M, A(2) , ..., A(N) , 

C L(l), ..., L(M) respectively. (Note that the A's are read 

C in degrees and then converted into radians. Also note that, 

C for an isotropic flux, N=1 and no A's are listed in 
C HFLUX00.DAT.) The file HFLUX01.DAT contains the values of 
C HA (1,1) , HA (1,2) , ..., HA ( 1 , M) where HA(J,K) is the integral 

C (in LET) directional flux ( in l/m**2-sec-sr) evaluated at 
C LET=L(K) and in the direction corresponding to the JTH 
C angular bin. Etc. for the other files. Note that N must 
C not exceed 30 and M must not exceed 200. After reading HA, 

C convert to the units l/cm**2— day— sr . 

PI=ATAN2 (0.0, -1.0) 

PI2=ATAN2 (1.0, 0.0) 

WRITE (*,*)' ENTER NORMAL INCIDENT CROSS SECTION (CM**2) 
READ*, AC 
R=SQRT (AC /PI ) 

WRITE (*,*) 'ENTER THRESHOLD LET (MEV-CM**2/MG) ' 

READ*, LTH 

WRITE (*,*) 'ENTER THICKNESS (CM)' 

READ* , T 

OPEN (UNIT=8 , STATUS= ' OLD ' , FILE= ' HFLUX00 . DAT ' ) 

REWIND (8) 

READ ( 8 , * ) N 
READ ( 8 , * ) M 
A ( 1) =0 . 0 
A (N+l) =PI 
IF (N.GE.2) THEN 
DO 10 J=2 , N 
READ ( 8 , * ) A ( J) 

A ( J) =A ( J) *PI/ 180.0 
10 CONTINUE 

END IF 
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PRECEDING PAGE BLANK NOT FILMED 



DO 20 J=1,M 
READ (8,*) L ( J) 

20 CONTINUE 
CLOSE (8) 

DO 4 0 J=1 , N 
JC1=M0D(J, 10) 

JC2= ( J-JC1) /10 

FLNM1=CHAR ( JC2+48 ) / /CHAR ( JC1+48 ) 

FLNM2= ' HFLUX ' / / FLNM1 / / ' . DAT ' 

OPEN (UNIT=8 , STATUS= • OLD 1 , FILE=FLNM2 ) 

REWIND (8) 

DO 30 K=1 , M 
READ ( 8 , * ) HA ( J, K) 

HA ( J, K) =8 . 64 *HA ( J , K) 

30 CONTINUE 
CLOSE (8) 

40 CONTINUE 

C Construct the B's using (13), the C's using (16), 

C and the D's using (22) . 

B ( 1 ) =0 . 0 

B ( 2 ) =ATAN ( 2 . 0*R/T) 

B ( 3 ) =PI2 

B (4 ) =PI-ATAN ( 2 . 0*R/T) 

B (5) =PI 
DO 60 1=1,4 
DO 50 J=1 , N 

Cl (I, J)=AMAX1(A(J) , B ( I ) ) 

C2 (I, J)=AMIN1(A(J+1) , B ( 1+1) ) 

50 CONTINUE 
60 CONTINUE 

DO 70 J=1 , N 

D1 ( 1 , J) =COS (Cl ( 1 , J) ) 

D1 (2 , J) =COS (Cl (2 , J) ) 

D1 ( 3 , J) =COS (PI-C1(3,J) ) 

D1 ( 4 , J) =COS (PI -Cl (4 , J) ) 

D2 (1 , J) =COS (C2 (1 , J) ) 

D2 (2 , J) =COS (C2 (2 , J) ) 

D2 (3 , J) =COS (PI-C2 ( 3 , J) ) 

D2 (4 , J) =COS (PI— C2 (4, J) ) 

70 CONTINUE 

C Construct S(K) using (32) and the F*s using (33). The W's 
C hold intermediate results and are reused in different 
C calculations. 

DO 80 K=1,M 

W=2 . 0*R*L (K) / (T*LTH) 

W=1.0-W*W 
S (K) =0 . 0 

IF (W.GT.0.0) S (K) =SQRT (W) 

F3 (K) =0 . 0 

IF (L (K) .LE.LTH) THEN 
W=L(K) /LTH 
W1=2.0*R*R*W*W 
W2=R*T*W*SQRT ( 1 . 0-W*W) 

W3=R*T*ASIN(W) 

W4= (4 . 0/3 . 0) *R*T* (LTH/L(K) ) *(1.0-W*W) **1.5 
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F3 (K)=W1+W2+W3+W4 
END IF 
W=S(K) 

Wl=2 . 0*R*R*W*W 
W2=R*T*W*SQRT ( 1 . 0-W*W) 

W3=R*T*ASIN(W) 

W4=(4 . 0/3 . 0) *R*T* ( LTH / L ( K ) ) * (1. 0-W*W) **1.5 
F4 (K) =W1+W2+W3+W4 
80 CONTINUE 

DO 110 1=1,4 
DO 100 J=1,N 
DO 90 K=1 , M 
W=D1 ( I , J) 

W1=2.0*R*R*W*W 
W2=R*T*W*SQRT ( 1 . 0-W*W) 

W3=R*T*ASIN(W) 

W4 = ( 4 . 0/ 3 . 0 ) *R*T* (LTH/L (K) ) * ( 1 . 0-W*W) **1 . 5 
FI ( I , J , K) =W1+W2+W3+W4 
W=D2 (I, J) 

W1=2.0*R*R*W*W 
W2=R*T*W*SQRT ( 1 . 0-W*W) 

W3=R*T*ASIN (W) 

W4 = (4 . 0/3 . 0) *R*T* (LTH/L (K) ) * ( 1 . 0-W*W) **1 . 5 
F2 (I, J,K)=W1+W2+W3+W4 
90 CONTINUE 
100 CONTINUE 
110 CONTINUE 

C Now calculate G(J,K) using (31). The W's hold 
C intermediate results. 

DO 130 J=1 , N 
DO 120 K=1,M 
W=L (K) /LTH 
U=S (K) 

W1=0. 0 

IF ( (C2 ( 1 , J) . GT . Cl ( 1 , J) ) . AND. (W. GT.D2 (1, J) ) ) W1=F3 (K) -F2 ( 1 , J, K) 
W2=0.0 

IF ( (C2 (1, J) .GT.C1(1, J) ) .AND. (W.GT.D1(1, J) ) ) W2=F1 ( 1 , J , K) -F3 (K) 
W3=0. 0 

IF ( (C2 ( 2 , J) . GT . Cl ( 2 , J) ) . AND . (D2 ( 2 , J) . GT . U) ) W3=F4 (K) -F2 (2 , J, K) 
W4=0 . 0 

IF ( (C2 (2 , J) . GT. Cl ( 2 , J) ) .AND . (D1 (2 , J) .GT. U) ) W4=F1 (2 , J,K) -F4 (K) 
W5=0 . 0 

IF ( (C2 ( 3 , J) . GT . Cl ( 3 , J) ) .AND . (D1 (3 , J) .GT.U) ) W5=F4 (K) -FI ( 3 , J , K) 
W6=0 . 0 

IF ( (C2(3,J) .GT.C1(3,J) ) .AND. (D2(3,J) .GT.U) ) W6=F2 ( 3 , J, K) -F4 (K) 
W7=0 . 0 

IF ( (C2 (4 , J) . GT . Cl (4 , J) ) .AND. (W . GT .D1(4,J) ) ) W7=F3 (K) -FI ( 4 , J, K) 
W8=0 . 0 

IF ( (C2 (4, J) .GT.C1(4, J) ) .AND. (W.GT.D2 (4, J) ) ) W8=F2 (4 , J, K) -F3 (K) 
G ( J, K) =W1+W2+W3+W4+W5+W6+W7+W8 
120 CONTINUE 
130 CONTINUE 

C Calculate the sum in (29) and print the result. 

RATE=0 . 0 
DO 150 J=1 , N 
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DO 140 K=1 , M-l 

RATE=RATE+ (G ( J, K+l) +G ( J, K) ) * (HA ( J , K) -HA ( J, K+l ) ) 
140 CONTINUE 
150 CONTINUE 

RATE=PI*RATE 
WRITE ( * , * ) ' LTH = ' , LTH 
WRITE (*, *) 'T = ' ,T 
WRITE (*, *) 'AC = ' ,AC 

WRITE (*,*) 'RATE = ' , RATE, ' UPSETS/DAY' 

END 
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APPENDIX B 


PROGRAM EFFLUX 

C Copyright (c) 1991, California Institute of Technology. 

C U.S. Government Sponsorship under NASA Contract NAS7-918 
C is acknowledged. 

C 

C This program is the same as NIFUR except that threshold 
C LET (LTH) , normal incident cross section (AC) , and 
C thickness (T) are not assigned by the user. AC is set 
C equal to 1 cm**2 and LTH is assigned in a loop so that it 
C is given each of the LET values that the flux is tabulated 
C against. For each LTH, T is varied for a maximum upset rate. 
C A tabulation of maximum rate versus LTH is the effective 
C flux and is stored in the file EFFLUX. OUT. 

REAL LTH, L( 200) 

CHARACTER* 2 , FLNM1 
CHARACTER* 1 1 , FLNM2 

DIMENSION A ( 3 1) , B ( 5 ) , Cl (4, 30) ,C2(4,30),D1(4,30) ,D2(4,30) 
DIMENSION S (200) ,RAT(200) ,TT(200) ,HA(30,200) 

DIMENSION F3 (200) , F4 (200) ,F1(4,30,200) ,F2(4,30,200) 
DIMENSION G (30 , 200) 

C Read data. The files HFLUX00.DAT, HFLUX01.DAT, ..., 
CHFLUX09.DAT, HFLUX10.DAT, HFLUX11.DAT, ... contain 
C environmental data. There is one file (after the "00" file) 

C for each angular bin (discussed below) . The angular inter- 
C val [0,180 degrees) is partitioned into the bins 
C [A(l) ,A(2) ) , [ A ( 2 ) ,A(3) ) , . . . , [ A (N) , A (N+l) ) where A(1)=0 

C and A (N+l) =180 . The discrete values L(l), ...,L(M) of LET 

C (in MeV-cm**2/mg) are selected by the user. The file 
C HFLUX00.DAT contains the values of N, M, A (2) , ..., A(N) , 

C L(l) , ..., L(M) respectively. (Note that the A's are read 
C in degrees and then converted into radians. Also note that, 

C for an isotropic flux, N=1 and no A's are listed in 
C HFLUX00.DAT.) The file HFLUX01.DAT contains the values of 
C HA (1,1) , HA (1,2) , ..., HA ( 1 , M) where HA(J,K) is the integral 

C (in LET) directional flux ( in l/m**2-sec-sr) evaluated at 
C LET=L (K) and in the direction corresponding to the JTH 
C angular bin. Etc. for the other files. Note that N must 
C not exceed 30 and M must not exceed 200. After reading HA, 

C convert to the units l/cm**2-day-sr . 

PI=ATAN2 ( 0 . 0 , -1 . 0) 

PI2=ATAN2 ( 1 . 0 , 0 . 0 ) 

AC=1 . 0 

R=SQRT (AC/PI) 

OPEN (UNIT=8 , STATUS= ' OLD ' , FILE= ' HFLUX00 . DAT ' ) 

REWIND ( 8 ) 

READ ( 8 , * ) N 
READ ( 8 , * ) M 
A ( 1) =0 . 0 
A (N+l) =PI 
IF (N.GE.2) THEN 
DO 10 J=2 , N 
READ (8,*) A ( J) 

A ( J) =A ( J) *PI / 180 . 0 


49 



10 CONTINUE 

END IF 
DO 20 J=1 , M 
READ (8,*) L ( J) 

20 CONTINUE 
CLOSE (8) 

DO 40 J=1 , N 
JC1=M0D(J, 10) 

JC2=(J-JC1) /10 

FLNM1=CHAR ( JC2+48 ) / /CHAR ( JC1+48 ) 

FLNM2 = ' HFLUX ' / / FLNM1 / / ' . DAT ' 

OPEN (UNIT=8 , STATUS= ' OLD ' , FILE=FLNM2 ) 

REWIND (8) 

DO 30 K=1 , M 
READ ( 8 , * ) HA ( J, K) 

HA ( J, K) =8 . 64*HA ( J, K) 

30 CONTINUE 
CLOSE (8) 

40 CONTINUE 

C Start an outer loop in J1 that assigns values to LTH. 

C An inner loop in J2 assigns values to T and calulates the 
C upset rate. The maximum rate calculated in the inner loop 
C is stored in RAT(Jl) and the value of T producing this 
C rate is stored in TT(Jl) . 

DO 210 Jl=l , M 
LTH=L( Jl) 

TT ( Jl) =0 . 0 
RAT ( Jl) =0 . 0 
DO 200 J2=l , 2 1 

T= 1 0 . 0 * * ( ( FLOAT (J2)-21.0)/10.0) 

C Now the upset rate calculation starts. 

C Construct the B's using (13), the C's using (16), 

C and the D's using (22). 

B ( 1) =0 . 0 

B ( 2 ) =ATAN ( 2 . 0*R/T) 

B ( 3 ) =PI2 

B (4 ) =PI-ATAN (2 . 0*R/T) 

B ( 5 ) =PI 
DO 60 1=1,4 
DO 50 J=1 , N 

Cl (I, J)=AMAX1(A(J) , B ( I ) ) 

C2 ( I , J) =AMIN1 (A ( J+l) , B ( 1+1) ) 

50 CONTINUE 
60 CONTINUE 

DO 70 J=1 , N 

D1 ( 1 , J) =COS (Cl ( 1 , J) ) 

D1 (2 , J) =COS (Cl (2 , J) ) 

D1 ( 3 , J) =COS (PI-C1 ( 3 , J) ) 

D1 (4 , J) =COS (PI -Cl (4, J) ) 

D2 ( 1 , J) =COS (C2 ( 1 , J) ) 

D2 (2 , J) =COS (C2 (2 , J) ) 

D2 ( 3 , J) =COS (PI-C2 (3 , J) ) 

D2 ( 4 , J) =COS (PI-C2 ( 4 , J) ) 

70 CONTINUE 

C Construct S(K) using (32) and the F's using (33). The W's 
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C hold intermediate results and are reused in different 
C calculations. 

DO 80 K=1,M 

W=2 . 0*R*L (K) / (T*LTH) 

W=l. 0-W*W 
S (K) =0 . 0 

IF (W.GT.0.0) S (K) =SQRT (W) 

F3 (K) =0 . 0 

IF (L(K) .LE.LTH) THEN 
W=L (K) /LTH 
Wl=2 . 0*R*R*W*W 
W2=R*T*W*SQRT ( 1 . 0-W*W) 

W3=R*T*ASIN(W) 

W4= (4 . 0/3.0) *R*T* (LTH/L (K) ) * (1. 0-W*W) **1.5 
F3 (K) =W1+W2+W3+W4 
END IF 
W=S(K) 

Wl=2 . 0*R*R*W*W 
W2=R*T*W*SQRT ( 1 . 0-W*W) 

W3=R*T*ASIN(W) 

W4= (4 . 0/3 . 0) *R*T* (LTH/L (K) ) * ( 1 . 0-W*W) **1 . 5 
F4 (K) =W1+W2+W3+W4 
80 CONTINUE 

DO 110 1=1,4 
DO 100 J=1 , N 
DO 90 K=1 , M 
W=D1 ( I , J) 

Wl=2 . 0*R*R*W*W 
W2=R*T*W*SQRT ( 1 . 0-W*W) 

W3=R*T*ASIN(W) 

W4=(4 . 0/3 . 0) *R*T* (LTH/L (K) ) *(1. 0-W*W) **1.5 
FI ( I , J , K) =W1+W2+W3+W4 
W=D2 (I, J) 

Wl=2 . 0*R*R*W*W 
W2=R*T*W*SQRT ( 1 . 0-W*W) 

W3=R*T*ASIN (W) 

W4= (4 .0/3.0) *R*T* (LTH/L (K) ) * (1 . 0-W*W) **1. 5 
F2 (I, J,K)=W1+W2+W3+W4 
90 CONTINUE 
100 CONTINUE 
110 CONTINUE 

C Now calculate G(J,K) using (31). The W's hold 
C intermediate results. 

DO 130 J=1 , N 
DO 120 K=1,M 
W=L (K) /LTH 
U=S(K) 

W1=0.0 

IF ( (C2 (1, J) .GT.C1(1, J) ) .AND. (W.GT.D2 (1, J) ) ) W1=F3 (K) -F2 ( 1 , J , K) 
W2=0 . 0 

IF ( (C2(1,J) .GT.C1(1,J) ) .AND. (W.GT.D1 (1, J) ) ) W2=F1 ( 1 , J , K) -F3 (K) 
W3=0. 0 

IF ( (C2 (2 , J) .GT. Cl (2 , J) ) .AND. (D2 (2, J) .GT.U) ) W3=F4 (K) -F2 (2 , J, K) 
W4=0 . 0 

IF ( (C2 (2, J) .GT.C1(2, J) ) .AND. (Dl(2, J) .GT.U) ) W4=F1 (2 , J,K) -F4 (K) 
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W5=0 . 0 

IF ( (C2 ( 3 , J) . GT . Cl ( 3 , J) ) .AND . (D1 ( 3 , J) . GT.U) ) W5=F4 (K) -FI ( 3 , J , K) 
W6=0 . 0 

IF ((C2(3,J) .GT.C1(3,J)) .AND. (D2(3,J) .GT.U)) W6=F2 (3,J,K) -F4 (K) 
W7=0 . 0 

IF ( (C2 (4 , J) . GT . Cl ( 4 , J) ) .AND. (W . GT . D1 (4 , J) ) ) W7=F3 (K) -FI (4 , J , K) 
W8=0 . 0 

IF ( (C2(4,J) ,GT.C1(4,J) ) .AND. (W . GT . D2 (4 , J) ) ) W8=F2 ( 4 , J, K) -F3 (K) 
G ( J , K) =W1+W2 +W3 +W4 +W5+W6 +W7 +W8 
120 CONTINUE 
130 CONTINUE 

C Calculate the upset rate using (29) . 

RATE=0 . 0 
DO 150 J=1 , N 
DO 140 K=1 , M-l 

RATE=RATE+ (G ( J, K+l) +G ( J, K) ) * (HA ( J , K) -HA( J, K+l) ) 

140 CONTINUE 
150 CONTINUE 

RATE=PI *RATE 

C Store the maximum upset rate and associated value of T 
C in RAT(Jl) and TT(J1). 

IF (RATE. GT. RAT ( Jl) ) THEN 
RAT(Jl) =RATE 
TT ( J 1 ) =T 
END IF 

200 CONTINUE 
210 CONTINUE 

C Record the data. The values of TT(J1) are of academic 
C interest and can be obtained by inserting additional 
C WRITE statements. 

OPEN (UNIT=8,STATUS=' UNKNOWN' , FI LE=' EFFLUX. OUT' ) 

REWIND (8) 

WRITE (8,*) 'EFFECTIVE FLUX' 

WRITE ( 8 , * ) ' ' 

WRITE (8, *) 'LET IN UNITS OF MeV-cm**2/mg ' 

WRITE(8,*) 'FLUX IN UNITS OF 1/ (cm**2-day) • 

WRITE (8, *) ' ' 

WRITE (8,230) ' LET ' , ' FLUX ' 

WRITE (8,230) ' ' , ' ' 

DO 220 J=1,M 

WRITE (8,240) L(J),RAT(J) 

220 CONTINUE 
CLOSE (8) 

230 FORMAT (A2 8, A20) 

240 FORMAT ( E3 1 . 4 , E2 0 . 4 ) 

END 
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